1//////////////////////////////////////////////////////////////////////
   2// LibFile: regions.scad
   3//   This file provides 2D Boolean set operations on polygons, where you can
   4//   compute, for example, the intersection or union of the shape defined by point lists, producing
   5//   a new point list.  Of course, such operations may produce shapes with multiple
   6//   components.  To handle that, we use "regions" which are lists of paths representing the polygons.
   7//   In addition to set operations, you can calculate offsets, determine whether a point is in a
   8//   region and you can decompose a region into parts.  
   9// Includes:
  10//   include <BOSL2/std.scad>
  11// FileGroup: Advanced Modeling
  12// FileSummary: Offsets and Boolean geometry of 2D paths and regions.
  13// FileFootnotes: STD=Included in std.scad
  14//////////////////////////////////////////////////////////////////////
  15
  16
  17// CommonCode:
  18//   include <BOSL2/rounding.scad>
  19
  20
  21// Section: Regions
  22//   A region is a list of polygons meeting these conditions:
  23//   .
  24//   - Every polygon on the list is simple, meaning it does not intersect itself
  25//   - Two polygons on the list do not cross each other
  26//   - A vertex of one polygon never meets the edge of another one except at a vertex
  27//   .
  28//   Note that this means vertex-vertex touching between two polygons is acceptable
  29//   to define a region.  Note, however, that regions with vertex-vertex contact usually
  30//   cannot be rendered with CGAL.  See {{is_valid_region()}} for examples of valid regions and
  31//   lists of polygons that are not regions.  Note that {{is_region_simple()}} will identify
  32//   regions with no polygon intersections at all, which should render successfully witih CGAL.  
  33//   .
  34//   The actual geometry of the region is defined by XORing together
  35//   all of the polygons in the list.  This may sound obscure, but it simply means that nested
  36//   boundaries make rings in the obvious fashion, and non-nested shapes simply union together.
  37//   Checking that a list of polygons is a valid region, meaning that it satisfies all of the conditions
  38//   above, can be a time consuming test, so it is not done automatically.  It is your responsibility to ensure that your regions are
  39//   compliant.  You can construct regions by making a suitable list of polygons, or by using
  40//   set operation function such as union() or difference(), which all acccept polygons, as
  41//   well as regions, as their inputs.  And if you must you can clean up an ill-formed region using make_region(),
  42//   which will break up self-intersecting polygons and polygons that cross each other.  
  43
  44
  45// Function: is_region()
  46// Synopsis: Returns true if the input appears to be a region.
  47// Topics: Regions, Paths, Polygons, List Handling
  48// See Also: is_valid_region(), is_1region(), is_region_simple()
  49// Usage:
  50//   bool = is_region(x);
  51// Description:
  52//   Returns true if the given item looks like a region.  A region is a list of non-crossing simple polygons.  This test just checks
  53//   that the argument is a list whose first entry is a path.  
  54function is_region(x) = is_list(x) && is_path(x.x);
  55
  56
  57// Function: is_valid_region()
  58// Synopsis: Returns true if the input is a valid region.
  59// Topics: Regions, Paths, Polygons, List Handling
  60// See Also: is_region(), is_1region(), is_region_simple()
  61// Usage:
  62//   bool = is_valid_region(region, [eps]);
  63// Description:
  64//   Returns true if the input is a valid region, meaning that it is a list of simple polygons whose segments do not cross each other.
  65//   This test can be time consuming with regions that contain many points.
  66//   It differs from `is_region()` which simply checks that the object is a list whose first entry is a path
  67//   because it searches all the list polygons for any self-intersections or intersections with each other.  
  68//   Will also return true if given a single simple polygon.  Use {{make_region()}} to convert sets of self-intersecting polygons into
  69//   a region.  
  70// Arguments:
  71//   region = region to check
  72//   eps = tolerance for geometric comparisons.  Default: `EPSILON` = 1e-9
  73// Example(2D,NoAxes):  In all of the examples each polygon in the region appears in a different color.  Two non-intersecting squares make a valid region.
  74//   region = [square(10), right(11,square(8))];
  75//   rainbow(region)stroke($item, width=.2,closed=true);
  76//   back(11)text(is_valid_region(region) ? "region" : "non-region", size=2);
  77// Example(2D,NoAxes):  Nested squares form a region
  78//   region = [for(i=[3:2:10]) square(i,center=true)];
  79//   rainbow(region)stroke($item, width=.2,closed=true);
  80//   back(6)text(is_valid_region(region) ? "region" : "non-region", size=2,halign="center");
  81// Example(2D,NoAxes):  Also a region:
  82//   region= [square(10,center=true), square(5,center=true), right(10,square(7))];
  83//   rainbow(region)stroke($item, width=.2,closed=true);
  84//   back(8)text(is_valid_region(region) ? "region" : "non-region", size=2);
  85// Example(2D,NoAxes):  The squares cross each other, so not a region
  86//   object = [square(10), move([8,8], square(8))];
  87//   rainbow(object)stroke($item, width=.2,closed=true);
  88//   back(17)text(is_valid_region(object) ? "region" : "non-region", size=2);
  89// Example(2D,NoAxes): A union is one way to fix the above example and get a region.  (Note that union is run here on two simple polygons, which are valid regions themselves and hence acceptable inputs to union.
  90//   region = union([square(10), move([8,8], square(8))]);
  91//   rainbow(region)stroke($item, width=.25,closed=true);
  92//   back(12)text(is_valid_region(region) ? "region" : "non-region", size=2);
  93// Example(2D,NoAxes):  Not a region due to a self-intersecting (non-simple) hourglass polygon
  94//   object = [move([-2,-2],square(14)), [[0,0],[10,0],[0,10],[10,10]]];
  95//   rainbow(object)stroke($item, width=.2,closed=true);
  96//   move([-1.5,13])text(is_valid_region(object) ? "region" : "non-region", size=2);
  97// Example(2D,NoAxes):  Breaking hourglass in half fixes it.  Now it's a region:
  98//   region = [move([-2,-2],square(14)), [[0,0],[10,0],[5,5]], [[5,5],[0,10],[10,10]]];
  99//   rainbow(region)stroke($item, width=.2,closed=true);
 100// Example(2D,NoAxes):  A single polygon corner touches an edge, so not a region:
 101//   object = [[[-10,0], [-10,10], [20,10], [20,-20], [-10,-20],
 102//              [-10,-10], [0,0], [10,-10], [10,0]]];
 103//   rainbow(object)stroke($item, width=.3,closed=true);
 104//   move([-4,12])text(is_valid_region(object) ? "region" : "non-region", size=3);
 105// Example(2D,NoAxes):  Corners touch in the same polygon, so the polygon is not simple and the object is not a region.
 106//   object = [[[0,0],[10,0],[10,10],[-10,10],[-10,0],[0,0],[-5,5],[5,5]]];
 107//   rainbow(object)stroke($item, width=.3,closed=true);
 108//   move([-10,12])text(is_valid_region(object) ? "region" : "non-region", size=3);
 109// Example(2D,NoAxes):  The shape above as a valid region with two polygons:
 110//   region = [  [[0,0],[10,0],[10,10],[-10,10],[-10,0]],
 111//               [[0,0],[5,5],[-5,5]]  ];
 112//   rainbow(region)stroke($item, width=.3,closed=true);
 113//   move([-5.5,12])text(is_valid_region(region) ? "region" : "non-region", size=3);
 114// Example(2D,NoAxes):  As with the "broken" hourglass, Touching at corners is OK.  This is a region.
 115//   region = [square(10), move([10,10], square(8))];
 116//   rainbow(region)stroke($item, width=.25,closed=true);
 117//   back(12)text(is_valid_region(region) ? "region" : "non-region", size=2);
 118// Example(2D,NoAxes): These two squares share part of an edge, hence not a region
 119//   object = [square(10), move([10,2], square(7))];
 120//   stroke(object[0], width=0.2,closed=true);
 121//   color("red")dashed_stroke(object[1], width=0.25,closed=true);
 122//   back(12)text(is_valid_region(object) ? "region" : "non-region", size=2);
 123// Example(2D,NoAxes): These two squares share a full edge, hence not a region
 124//   object = [square(10), right(10, square(10))];
 125//   stroke(object[0], width=0.2,closed=true);
 126//   color("red")dashed_stroke(object[1], width=0.25,closed=true);
 127//   back(12)text(is_valid_region(object) ? "region" : "non-region", size=2);
 128// Example(2D,NoAxes): Sharing on edge on the inside, also not a regionn
 129//   object = [square(10), [[0,0], [2,2],[2,8],[0,10]]];
 130//   stroke(object[0], width=0.2,closed=true);
 131//   color("red")dashed_stroke(object[1], width=0.25,closed=true);
 132//   back(12)text(is_valid_region(object) ? "region" : "non-region", size=2);
 133// Example(2D,NoAxes): Crossing at vertices is also bad
 134//   object = [square(10), [[10,0],[0,10],[8,13],[13,8]]];
 135//   rainbow(object)stroke($item, width=.2,closed=true);
 136//   back(14)text(is_valid_region(object) ? "region" : "non-region", size=2);
 137// Example(2D,NoAxes): One polygon touches another in the middle of an edge
 138//   object = [square(10), [[10,5],[15,0],[15,10]]];
 139//   rainbow(object)stroke($item, width=.2,closed=true);
 140//   back(11)text(is_valid_region(object) ? "region" : "non-region", size=2);
 141// Example(2D,NoAxes): The polygon touches the side, but the side has a vertex at the contact point so this is a region
 142//   poly1 = [ each square(30,center=true), [15,0]];
 143//   poly2 = right(10,circle(5,$fn=4));
 144//   poly3 = left(0,circle(5,$fn=4));
 145//   poly4 = move([0,-8],square([10,3]));
 146//   region = [poly1,poly2,poly3,poly4];
 147//   rainbow(region)stroke($item, width=.25,closed=true);
 148//   move([-5,16.5])text(is_valid_region(region) ? "region" : "non-region", size=3);
 149//   color("black")move_copies(region[0]) circle(r=.4);
 150// Example(2D,NoAxes): The polygon touches the side, but not at a vertex so this is not a region
 151//   poly1 = fwd(4,[ each square(30,center=true), [15,0]]);
 152//   poly2 = right(10,circle(5,$fn=4));
 153//   poly3 = left(0,circle(5,$fn=4));
 154//   poly4 = move([0,-8],square([10,3]));
 155//   object = [poly1,poly2,poly3,poly4];
 156//   rainbow(object)stroke($item, width=.25,closed=true);
 157//   move([-9,12.5])text(is_valid_region(object) ? "region" : "non-region", size=3);
 158//   color("black")move_copies(object[0]) circle(r=.4);
 159// Example(2D,NoAxes): The inner polygon touches the middle of the edges, so not a region
 160//   poly1 = square(20,center=true);
 161//   poly2 = circle(10,$fn=8);
 162//   object=[poly1,poly2];
 163//   rainbow(object)stroke($item, width=.25,closed=true);
 164//   move([-10,11.4])text(is_valid_region(object) ? "region" : "non-region", size=3);
 165// Example(2D,NoAxes): The above shape made into a region using {{difference()}} now has four components that touch at corners
 166//   poly1 = square(20,center=true);
 167//   poly2 = circle(10,$fn=8);
 168//   region = difference(poly1,poly2);
 169//   rainbow(region)stroke($item, width=.25,closed=true);
 170//   move([-5,11.4])text(is_valid_region(region) ? "region" : "non-region", size=3);
 171function is_valid_region(region, eps=EPSILON) =
 172   let(region=force_region(region))
 173   assert(is_region(region), "Input is not a region")
 174   // no short paths
 175   [for(p=region) if (len(p)<3) 1] == []
 176   &&
 177   // all paths are simple
 178   [for(p=region) if (!is_path_simple(p,closed=true,eps=eps)) 1] == []
 179   &&
 180   // paths do not cross each other
 181   [for(i=[0:1:len(region)-2])
 182            if (_polygon_crosses_region(list_tail(region,i+1),region[i], eps=eps)) 1] == []
 183   &&
 184   // one path doesn't touch another in the middle of an edge
 185   [for(i=idx(region), j=idx(region))
 186       if (i!=j) for(v=region[i], edge=pair(region[j],wrap=true))
 187           let(
 188               v1 = edge[1]-edge[0],
 189               v0 = v - edge[0],
 190               t = v0*v1/(v1*v1)
 191           )
 192           if (abs(cross(v0,v1))<eps*norm(v1) && t>eps && t<1-eps) 1
 193   ]==[];
 194
 195
 196
 197// internal function:
 198// returns true if the polygon crosses the region so that part of the 
 199// polygon is inside the region and part is outside.  
 200function _polygon_crosses_region(region, poly, eps=EPSILON) =
 201    let(  
 202        subpaths = flatten(split_region_at_region_crossings(region,[poly],eps=eps)[1])
 203    )
 204    [for(path=subpaths)
 205      let(isect=
 206         [for (subpath = subpaths)
 207          let(
 208                midpt = mean([subpath[0], subpath[1]]),
 209                rel = point_in_region(midpt,region,eps=eps)
 210          )
 211          rel
 212         ])
 213       if (!all_equal(isect) || isect[0]==0) 1 ] != [];
 214
 215
 216// Function: is_region_simple()
 217// Synopsis: Returns true if the input is a region with no corner contact.
 218// Topics: Regions, Paths, Polygons, List Handling
 219// See Also: is_region(), is_valid_region(), is_1region()
 220// Usage:
 221//   bool = is_region_simple(region, [eps]);
 222// Description:
 223//   We extend the notion of the simple path to regions: a simple region is entirely
 224//   non-self-intersecting, meaning that it is formed from a list of simple polygons that
 225//   don't intersect each other at all&mdash;not even with corner contact points.
 226//   Regions with corner contact are valid but may fail CGAL.  Simple regions
 227//   should not create problems with CGAL.  
 228// Arguments:
 229//   region = region to check
 230//   eps = tolerance for geometric comparisons.  Default: `EPSILON` = 1e-9
 231// Example(2D,NoAxes):  Corner contact means it's not simple
 232//   region = [move([-2,-2],square(14)), [[0,0],[10,0],[5,5]], [[5,5],[0,10],[10,10]]];
 233//   rainbow(region)stroke($item, width=.2,closed=true);
 234//   move([-1,13])text(is_region_simple(region) ? "simple" : "not-simple", size=2);
 235// Example(2D,NoAxes):  Moving apart the triangles makes it simple:
 236//   region = [move([-2,-2],square(14)), [[0,0],[10,0],[5,4.5]], [[5,5.5],[0,10],[10,10]]];
 237//   rainbow(region)stroke($item, width=.2,closed=true);
 238//   move([1,13])text(is_region_simple(region) ? "simple" : "not-simple", size=2);
 239function is_region_simple(region, eps=EPSILON) =
 240   let(region=force_region(region))
 241   assert(is_region(region), "Input is not a region")
 242   [for(p=region) if (!is_path_simple(p,closed=true,eps=eps)) 1] == []
 243   &&
 244   [for(i=[0:1:len(region)-2])
 245       if (_region_region_intersections([region[i]], list_tail(region,i+1), eps=eps)[0][0] != []) 1
 246   ] ==[];
 247  
 248  
 249// Function: make_region()
 250// Synopsis: Converts lists of intersecting polygons into valid regions.
 251// SynTags: Region
 252// Topics: Regions, Paths, Polygons, List Handling
 253// See Also: force_region(), region()
 254// 
 255// Usage:
 256//   region = make_region(polys, [nonzero], [eps]);
 257// Description:
 258//   Takes a list of polygons that may intersect themselves or cross each other 
 259//   and converts it into a properly defined region without these defects.
 260// Arguments:
 261//   polys = list of polygons to use
 262//   nonzero = set to true to use nonzero rule for polygon membership.  Default: false
 263//   eps = Epsilon for geometric comparisons.  Default: `EPSILON` (1e-9)
 264// Example(2D,NoAxes):  The pentagram is self-intersecting, so it is not a region.  Here it becomes five triangles:
 265//   pentagram = turtle(["move",100,"left",144], repeat=4);
 266//   region = make_region(pentagram);
 267//   rainbow(region)stroke($item, width=1,closed=true);
 268// Example(2D,NoAxes):  Alternatively with the nonzero option you can get the perimeter:
 269//   pentagram = turtle(["move",100,"left",144], repeat=4);
 270//   region = make_region(pentagram,nonzero=true);
 271//   rainbow(region)stroke($item, width=1,closed=true);
 272// Example(2D,NoAxes):  Two crossing squares become two L-shaped components
 273//   region = make_region([square(10), move([5,5],square(8))]);
 274//   rainbow(region)stroke($item, width=.3,closed=true);
 275
 276function make_region(polys,nonzero=false,eps=EPSILON) =
 277     let(polys=force_region(polys))
 278     assert(is_region(polys), "Input is not a region")
 279     exclusive_or(
 280                  [for(poly=polys) each polygon_parts(poly,nonzero,eps)],
 281                  eps=eps);
 282
 283// Function: force_region()
 284// Synopsis: Given a polygon returns a region.
 285// SynTags: Region
 286// Topics: Regions, Paths, Polygons, List Handling
 287// See Also: make_region(), region()
 288// Usage:
 289//   region = force_region(poly)
 290// Description:
 291//   If the input is a polygon then return it as a region.  Otherwise return it unaltered.
 292// Arguments:
 293//   poly = polygon to turn into a region
 294function force_region(poly) = is_path(poly) ? [poly] : poly;
 295
 296
 297// Section: Turning a region into geometry
 298
 299// Module: region()
 300// Synopsis: Creates the 2D polygons described by the given region or list of polygons.
 301// SynTags: Geom
 302// Topics: Regions, Paths, Polygons, List Handling
 303// See Also: make_region(), debug_region()
 304// Usage:
 305//   region(r, [anchor], [spin=], [cp=], [atype=]) [ATTACHMENTS];
 306// Description:
 307//   Creates the 2D polygons described by the given region or list of polygons.  This module works on
 308//   arbitrary lists of polygons that cross each other and hence do not define a valid region.  The
 309//   displayed result is the exclusive-or of the polygons listed in the input. 
 310// Arguments:
 311//   r = region to create as geometry
 312//   anchor = Translate so anchor point is at origin (0,0,0).  See [anchor](attachments.scad#subsection-anchor).  Default: `"origin"`
 313//   ---
 314//   spin = Rotate this many degrees after anchor.  See [spin](attachments.scad#subsection-spin).  Default: `0`
 315//   cp = Centerpoint for determining intersection anchors or centering the shape.  Determintes the base of the anchor vector.  Can be "centroid", "mean", "box" or a 2D point.  Default: "centroid"
 316//   atype = Set to "hull" or "intersect" to select anchor type.  Default: "hull"
 317// Example(2D): Displaying a region
 318//   region([circle(d=50), square(25,center=true)]);
 319// Example(2D): Displaying a list of polygons that intersect each other, which is not a region
 320//   rgn = concat(
 321//       [for (d=[50:-10:10]) circle(d=d-5)],
 322//       [square([60,10], center=true)]
 323//   );
 324//   region(rgn);
 325module region(r, anchor="origin", spin=0, cp="centroid", atype="hull")
 326{
 327    assert(in_list(atype, _ANCHOR_TYPES), "Anchor type must be \"hull\" or \"intersect\"");
 328    r = force_region(r);
 329    dummy=assert(is_region(r), "Input is not a region");
 330    points = flatten(r);
 331    lengths = [for(path=r) len(path)];
 332    starts = [0,each cumsum(lengths)];
 333    paths = [for(i=idx(r)) count(s=starts[i], n=lengths[i])];
 334    attachable(anchor, spin, two_d=true, region=r, extent=atype=="hull", cp=cp){
 335      polygon(points=points, paths=paths);
 336      children();
 337    }
 338}
 339
 340
 341
 342// Module: debug_region()
 343// Synopsis: Draws an annotated region.
 344// SynTags: Geom
 345// Topics: Shapes (2D)
 346// See Also: region(), debug_polygon(), debug_vnf(), debug_bezier()
 347//
 348// Usage:
 349//   debug_region(region, [vertices=], [edges=], [convexity=], [size=]);
 350// Description:
 351//   A replacement for {{region()}} that displays the region and labels the vertices and
 352//   edges.  The region vertices and edges are labeled with letters to identify the path
 353//   component in the region, starting with A.  
 354//   The start of each path is marked with a blue circle and the end with a pink diamond.
 355//   You can suppress the display of vertex or edge labeling using the `vertices` and `edges` arguments.
 356// Arguments:
 357//   region = region to display
 358//   ---
 359//   vertices = if true display vertex labels and start/end markers.  Default: true
 360//   edges = if true display edge labels.  Default: true
 361//   convexity = The max number of walls a ray can pass through the given polygon paths.
 362//   size = The base size of the line and labels.
 363// Example(2D,Big):
 364//   region = make_region([square(15), move([5,5],square(15))]);
 365//   debug_region(region,size=1);
 366module debug_region(region, vertices=true, edges=true, convexity=2, size=1)
 367{
 368  
 369  if (is_path(region) || (is_region(region) && len(region)==1))
 370    debug_polygon(force_path(region), vertices=vertices, edges=edges, convexity=convexity, size=size);
 371  else {
 372    for(i=idx(region))
 373      echo(str("points_",chr(97+i)," = ",region[i]))
 374    linear_extrude(height=0.01, convexity=convexity, center=true) 
 375      region(region);
 376    if(vertices)
 377        _debug_poly_verts(region,size);
 378    for(j=idx(region)){
 379      if(edges)
 380        _debug_poly_edges(j,region[j],vertices=vertices,size=size);
 381    }      
 382  }      
 383}
 384
 385
 386
 387// Section: Geometrical calculations with regions
 388
 389// Function: point_in_region()
 390// Synopsis: Tests if a point is inside, outside, or on the border of a region. 
 391// Topics: Regions, Points, Comparison
 392// See Also: region_area(), are_regions_equal()
 393// Usage:
 394//   check = point_in_region(point, region, [eps]);
 395// Description:
 396//   Tests if a point is inside, outside, or on the border of a region.  
 397//   Returns -1 if the point is outside the region.
 398//   Returns 0 if the point is on the boundary.
 399//   Returns 1 if the point lies inside the region.
 400// Arguments:
 401//   point = The point to test.
 402//   region = The region to test against, as a list of polygon paths.
 403//   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 404// Example(2D,Med):  Red points are in the region.
 405//   region = [for(i=[2:4:10]) hexagon(r=i)];
 406//   color("#ff7") region(region);
 407//   for(x=[-10:10], y=[-10:10])
 408//     if (point_in_region([x,y], region)>=0)
 409//       move([x,y]) color("red") circle(0.15, $fn=12);
 410//     else
 411//       move([x,y]) color("#ddf") circle(0.1, $fn=12);
 412function point_in_region(point, region, eps=EPSILON) =
 413    let(region=force_region(region))
 414    assert(is_region(region), "Region given to point_in_region is not a region")
 415    assert(is_vector(point,2), "Point must be a 2D point in point_in_region")
 416    _point_in_region(point, region, eps);
 417
 418function _point_in_region(point, region, eps=EPSILON, i=0, cnt=0) =
 419      i >= len(region) ? ((cnt%2==1)? 1 : -1)
 420    : let(
 421          pip = point_in_polygon(point, region[i], eps=eps)
 422      )
 423      pip == 0 ? 0
 424   : _point_in_region(point, region, eps=eps, i=i+1, cnt = cnt + (pip>0? 1 : 0));
 425
 426
 427// Function: region_area()
 428// Synopsis: Computes the area of the specified valid region.
 429// Topics: Regions, Area
 430// Usage:
 431//   area = region_area(region);
 432// Description:
 433//   Computes the area of the specified valid region. (If the region is invalid and has self intersections
 434//   the result is meaningless.)
 435// Arguments:
 436//   region = region whose area to compute
 437// Examples:
 438//   area = region_area([square(10), right(20,square(8))]);  // Returns 164
 439function region_area(region) =
 440  assert(is_region(region), "Input must be a region")
 441  let(
 442      parts = region_parts(region)
 443  )
 444  -sum([for(R=parts, poly=R) polygon_area(poly,signed=true)]);
 445
 446
 447
 448function _clockwise_region(r) = [for(p=r) clockwise_polygon(p)];
 449
 450// Function: are_regions_equal()
 451// Synopsis: Returns true if given regions are the same polygons.
 452// Topics: Regions, Polygons, Comparison
 453// Usage:
 454//    b = are_regions_equal(region1, region2, [either_winding])
 455// Description:
 456//    Returns true if the components of region1 and region2 are the same polygons (in any order). 
 457// Arguments:
 458//    region1 = first region
 459//    region2 = second region
 460//    either_winding = if true then two shapes test equal if they wind in opposite directions.  Default: false
 461function are_regions_equal(region1, region2, either_winding=false) =
 462    let(
 463        region1=force_region(region1),
 464        region2=force_region(region2)
 465    )
 466    assert(is_region(region1) && is_region(region2), "One of the inputs is not a region")
 467    len(region1) != len(region2)? false :
 468    __are_regions_equal(either_winding?_clockwise_region(region1):region1,
 469                        either_winding?_clockwise_region(region2):region2,
 470                        0);
 471
 472function __are_regions_equal(region1, region2, i) =
 473    i >= len(region1)? true :
 474    !_is_polygon_in_list(region1[i], region2)? false :
 475    __are_regions_equal(region1, region2, i+1);
 476
 477
 478/// Internal Function: _region_region_intersections()
 479/// Usage:
 480///    risect = _region_region_intersections(region1, region2, [closed1], [closed2], [eps]
 481/// Description:
 482///    Returns a pair of sorted lists such that risect[0] is a list of intersection
 483///    points for every path in region1, and similarly risect[1] is a list of intersection
 484///    points for the paths in region2.  For each path the intersection list is
 485///    a sorted list of the form [PATHIND, SEGMENT, U].  You can specify that the paths in either
 486///    region be regarded as open paths if desired.  Default is to treat them as
 487///    regions and hence the paths as closed polygons.
 488///    .
 489///    Included as intersection points are points where region1 touches itself at a vertex or
 490///    region2 touches itself at a vertex.  (The paths are assumed to have no self crossings.
 491///    Self crossings of the paths in the regions are not returned.)
 492function _region_region_intersections(region1, region2, closed1=true,closed2=true, eps=EPSILON) =
 493   let(
 494       intersections =   [
 495           for(p1=idx(region1))
 496              let(
 497                  path = closed1?list_wrap(region1[p1]):region1[p1]
 498              )
 499              for(i = [0:1:len(path)-2])
 500                  let(
 501                      a1 = path[i],
 502                      a2 = path[i+1],
 503                      nrm = norm(a1-a2)
 504                  )
 505                  if( nrm>eps )  // ignore zero-length path edges
 506                       let( 
 507                           seg_normal = [-(a2-a1).y, (a2-a1).x]/nrm,
 508                           ref = a1*seg_normal
 509                       )
 510                           // `signs[j]` is the sign of the signed distance from
 511                           // poly vertex j to the line [a1,a2] where near zero
 512                           // distances are snapped to zero;  poly edges 
 513                           //  with equal signs at its vertices cannot intersect
 514                           // the path edge [a1,a2] or they are collinear and 
 515                           // further tests can be discarded.
 516                       for(p2=idx(region2))
 517                           let(
 518                               poly  = closed2?list_wrap(region2[p2]):region2[p2],
 519                               signs = [for(v=poly*seg_normal) abs(v-ref) < eps ? 0 : sign(v-ref) ]
 520                           ) 
 521                           if(max(signs)>=0 && min(signs)<=0) // some edge intersects line [a1,a2]
 522                               for(j=[0:1:len(poly)-2]) 
 523                                   if(signs[j]!=signs[j+1])
 524                                        let( // exclude non-crossing and collinear segments
 525                                            b1 = poly[j],
 526                                            b2 = poly[j+1],
 527                                            isect = _general_line_intersection([a1,a2],[b1,b2],eps=eps) 
 528                                        )
 529                                        if (isect 
 530                                            && isect[1]>= -eps 
 531                                            && isect[1]<= 1+eps 
 532                                            && isect[2]>= -eps
 533                                            && isect[2]<= 1+eps)       
 534                                         [[p1,i,isect[1]], [p2,j,isect[2]]]
 535         ],
 536         regions=[region1,region2],
 537         // Create a flattened index list corresponding to the points in region1 and region2
 538         // that gives each point as an intersection point
 539         ptind = [for(i=[0:1])   
 540                    [for(p=idx(regions[i]))
 541                       for(j=idx(regions[i][p])) [p,j,0]]],
 542         points = [for(i=[0:1]) flatten(regions[i])],
 543         // Corner points are those points where the region touches itself, hence duplicate
 544         // points in the region's point set
 545         cornerpts = [for(i=[0:1])
 546                         [for(k=vector_search(points[i],eps,points[i]))
 547                             each if (len(k)>1) select(ptind[i],k)]],
 548         risect = [for(i=[0:1]) concat(column(intersections,i), cornerpts[i])],
 549         counts = [count(len(region1)), count(len(region2))],
 550         pathind = [for(i=[0:1]) search(counts[i], risect[i], 0)]
 551       )
 552       [for(i=[0:1]) [for(j=counts[i]) _sort_vectors(select(risect[i],pathind[i][j]))]];
 553         
 554
 555// Section: Breaking up regions into subregions
 556
 557
 558// Function: split_region_at_region_crossings()
 559// Synopsis: Splits regions where polygons touch and at intersections.
 560// Topics: Regions, Polygons, List Handling
 561// See Also: region_parts()
 562// 
 563// Usage:
 564//   split_region = split_region_at_region_crossings(region1, region2, [closed1], [closed2], [eps])
 565// Description:
 566//   Splits region1 at the places where polygons in region1 touches each other at corners and at locations
 567//   where region1 intersections region2.  Split region2 similarly with respect to region1.
 568//   The return is a pair of results of the form [split1, split2] where split1=[frags1,frags2,...]
 569//   and frags1 is a list of paths that when placed end to end (in the given order), give the first polygon of region1.
 570//   Each path in the list is either entirely inside or entirely outside region2.  
 571//   Then frags2 is the decomposition of the second polygon into path pieces, and so on.  Finally split2 is
 572//   the same list, but for the polygons in region2.  
 573//   You can pass a single polygon in for either region, but the output will be a singleton list, as if
 574//   you passed in a singleton region.  If you set the closed parameters to false then the region components
 575//   will be treated as open paths instead of polygons.  
 576// Arguments:
 577//   region1 = first region
 578//   region2 = second region
 579//   closed1 = if false then treat region1 as list of open paths.  Default: true
 580//   closed2 = if false then treat region2 as list of open paths.  Default: true
 581//   eps = Acceptable variance.  Default: `EPSILON` (1e-9)
 582// Example(2D): 
 583//   path = square(50,center=false);
 584//   region = [circle(d=80), circle(d=40)];
 585//   paths = split_region_at_region_crossings(path, region);
 586//   color("#aaa") region(region);
 587//   rainbow(paths[0][0]) stroke($item, width=2);
 588//   right(110){
 589//     color("#aaa") region([path]);
 590//     rainbow(flatten(paths[1])) stroke($item, width=2);
 591//   }
 592function split_region_at_region_crossings(region1, region2, closed1=true, closed2=true, eps=EPSILON) = 
 593    let(
 594        region1=force_region(region1),
 595        region2=force_region(region2)
 596    )
 597    assert(is_region(region1) && is_region(region2),"One of the inputs is not a region")
 598    let(
 599        xings = _region_region_intersections(region1, region2, closed1, closed2, eps),
 600        regions = [region1,region2],
 601        closed = [closed1,closed2]
 602    )
 603    [for(i=[0:1])
 604      [for(p=idx(xings[i]))
 605        let(
 606            crossings = deduplicate([
 607                                     [p,0,0],
 608                                     each xings[i][p],
 609                                     [p,len(regions[i][p])-(closed[i]?1:2), 1],
 610                                    ],eps=eps),
 611            subpaths = [
 612                for (frag = pair(crossings)) 
 613                    deduplicate(
 614                        _path_select(regions[i][p], frag[0][1], frag[0][2], frag[1][1], frag[1][2], closed=closed[i]),
 615                        eps=eps
 616                    )
 617            ]
 618        )
 619        [for(s=subpaths) if (len(s)>1) s]
 620       ]
 621    ];
 622                
 623                
 624
 625// Function: region_parts()
 626// Synopsis: Splits a region into a list of connected regions.
 627// SynTags: RegList
 628// Topics: Regions, List Handling
 629// See Also: split_region_at_region_crossings()
 630// Usage:
 631//   rgns = region_parts(region);
 632// Description:
 633//   Divides a region into a list of connected regions.  Each connected region has exactly one clockwise outside boundary
 634//   and zero or more counter-clockwise outlines defining internal holes.  Note that behavior is undefined on invalid regions whose
 635//   components cross each other.
 636// Example(2D,NoAxes):
 637//   R = [for(i=[1:7]) square(i,center=true)];
 638//   region_list = region_parts(R);
 639//   rainbow(region_list) region($item);
 640// Example(2D,NoAxes):
 641//   R = [back(7,square(3,center=true)),
 642//        square([20,10],center=true),
 643//        left(5,square(8,center=true)),
 644//        for(i=[4:2:8])
 645//          right(5,square(i,center=true))];
 646//   region_list = region_parts(R);
 647//   rainbow(region_list) region($item);
 648function region_parts(region) =
 649   let(
 650       region = force_region(region)
 651   )
 652   assert(is_region(region), "Input is not a region")
 653   let(
 654       inside = [for(i=idx(region))
 655                    let(pt = mean([region[i][0], region[i][1]]))
 656                    [for(j=idx(region))  i==j ? 0
 657                                       : point_in_polygon(pt,region[j]) >=0 ? 1 : 0]
 658                ],
 659       level = inside*repeat(1,len(region))
 660   )
 661   [ for(i=idx(region))
 662      if(level[i]%2==0)
 663         let(
 664             possible_children = search([level[i]+1],level,0)[0],
 665             keep=search([1], select(inside,possible_children), 0, i)[0]
 666         )
 667         [
 668           clockwise_polygon(region[i]),
 669           for(good=keep)
 670              ccw_polygon(region[possible_children[good]])
 671         ]
 672    ];
 673
 674
 675
 676
 677// Section: Offset and 2D Boolean Set Operations
 678
 679
 680function _offset_chamfer(center, points, delta) =
 681    is_undef(points[1])?
 682        let( points = select(points,[0,2]),
 683             center = mean(points),
 684             dir = sign(delta)*line_normal(points),
 685             halfside = tan(22.5)*abs(delta)
 686        )
 687        [ points[0]+dir*halfside,
 688          center + dir*abs(delta) + unit(points[0]-center)*halfside,
 689          center + dir*abs(delta) + unit(points[1]-center)*halfside, 
 690          points[1]+dir*halfside
 691        ]
 692    :
 693    let(
 694        dist = sign(delta)*norm(center-line_intersection(select(points,[0,2]), [center, points[1]])),
 695        endline = _shift_segment(select(points,[0,2]), delta-dist)
 696    ) [
 697        line_intersection(endline, select(points,[0,1])),
 698        line_intersection(endline, select(points,[1,2]))
 699    ];
 700
 701
 702function _shift_segment(segment, d) =
 703    assert(!approx(segment[0],segment[1]),"Path has repeated points")
 704    move(d*line_normal(segment),segment);
 705
 706
 707// Extend to segments to their intersection point.  First check if the segments already have a point in common,
 708// which can happen if two colinear segments are input to the path variant of `offset()`
 709function _segment_extension(s1,s2) =
 710    norm(s1[1]-s2[0])<1e-6 ? s1[1] : line_intersection(s1,s2,LINE,LINE);
 711
 712
 713function _makefaces(direction, startind, good, pointcount, closed) =
 714    let(
 715        lenlist = list_bset(good, pointcount),
 716        numfirst = len(lenlist),
 717        numsecond = sum(lenlist),
 718        prelim_faces = _makefaces_recurse(startind, startind+len(lenlist), numfirst, numsecond, lenlist, closed)
 719    )
 720    direction? [for(entry=prelim_faces) reverse(entry)] : prelim_faces;
 721
 722
 723function _makefaces_recurse(startind1, startind2, numfirst, numsecond, lenlist, closed, firstind=0, secondind=0, faces=[]) =
 724    // We are done if *both* firstind and secondind reach their max value, which is the last point if !closed or one past
 725    // the last point if closed (wrapping around).  If you don't check both you can leave a triangular gap in the output.
 726    ((firstind == numfirst - (closed?0:1)) && (secondind == numsecond - (closed?0:1)))? faces :
 727    _makefaces_recurse(
 728        startind1, startind2, numfirst, numsecond, lenlist, closed, firstind+1, secondind+lenlist[firstind],
 729        lenlist[firstind]==0? (
 730            // point in original path has been deleted in offset path, so it has no match.  We therefore
 731            // make a triangular face using the current point from the offset (second) path
 732            // (The current point in the second path can be equal to numsecond if firstind is the last point)
 733            concat(faces,[[secondind%numsecond+startind2, firstind+startind1, (firstind+1)%numfirst+startind1]])
 734            // in this case a point or points exist in the offset path corresponding to the original path
 735        ) : (
 736            concat(faces,
 737                // First generate triangular faces for all of the extra points (if there are any---loop may be empty)
 738                [for(i=[0:1:lenlist[firstind]-2]) [firstind+startind1, secondind+i+1+startind2, secondind+i+startind2]],
 739                // Finish (unconditionally) with a quadrilateral face
 740                [
 741                    [
 742                        firstind+startind1,
 743                        (firstind+1)%numfirst+startind1,
 744                        (secondind+lenlist[firstind])%numsecond+startind2,
 745                        (secondind+lenlist[firstind]-1)%numsecond+startind2
 746                    ]
 747                ]
 748            )
 749        )
 750    );
 751
 752
 753// Determine which of the shifted segments are good
 754function _good_segments(path, d, shiftsegs, closed, quality) =
 755    let(
 756        maxind = len(path)-(closed ? 1 : 2),
 757        pathseg = [for(i=[0:maxind]) select(path,i+1)-path[i]],
 758        pathseg_len =  [for(seg=pathseg) norm(seg)],
 759        pathseg_unit = [for(i=[0:maxind]) pathseg[i]/pathseg_len[i]],
 760        // Order matters because as soon as a valid point is found, the test stops
 761        // This order works better for circular paths because they succeed in the center
 762        alpha = concat([for(i=[1:1:quality]) i/(quality+1)],[0,1])
 763    ) [
 764        for (i=[0:len(shiftsegs)-1])
 765            (i>maxind)? true :
 766            _segment_good(path,pathseg_unit,pathseg_len, d - 1e-7, shiftsegs[i], alpha)
 767    ];
 768
 769
 770// Determine if a segment is good (approximately)
 771// Input is the path, the path segments normalized to unit length, the length of each path segment
 772// the distance threshold, the segment to test, and the locations on the segment to test (normalized to [0,1])
 773// The last parameter, index, gives the current alpha index.
 774//
 775// A segment is good if any part of it is farther than distance d from the path.  The test is expensive, so
 776// we want to quit as soon as we find a point with distance > d, hence the recursive code structure.
 777//
 778// This test is approximate because it only samples the points listed in alpha.  Listing more points
 779// will make the test more accurate, but slower.
 780function _segment_good(path,pathseg_unit,pathseg_len, d, seg,alpha ,index=0) =
 781    index == len(alpha) ? false :
 782    _point_dist(path,pathseg_unit,pathseg_len, alpha[index]*seg[0]+(1-alpha[index])*seg[1]) > d ? true :
 783    _segment_good(path,pathseg_unit,pathseg_len,d,seg,alpha,index+1);
 784
 785
 786// Input is the path, the path segments normalized to unit length, the length of each path segment
 787// and a test point.  Computes the (minimum) distance from the path to the point, taking into
 788// account that the minimal distance may be anywhere along a path segment, not just at the ends.
 789function _point_dist(path,pathseg_unit,pathseg_len,pt) =
 790    min([
 791        for(i=[0:len(pathseg_unit)-1]) let(
 792            v = pt-path[i],
 793            projection = v*pathseg_unit[i],
 794            segdist = projection < 0? norm(pt-path[i]) :
 795                projection > pathseg_len[i]? norm(pt-select(path,i+1)) :
 796                norm(v-projection*pathseg_unit[i])
 797        ) segdist
 798    ]);
 799
 800
 801// Function: offset()
 802// Synopsis: Takes a 2D path, polygon or region and returns a path offset by an amount.
 803// SynTags: Path, Region, Ext
 804// Topics: Paths, Polygons, Regions
 805// Usage:
 806//   offsetpath = offset(path, [r=|delta=], [chamfer=], [closed=], [check_valid=], [quality=], [same_length=])
 807//   path_faces = offset(path, return_faces=true, [r=|delta=], [chamfer=], [closed=], [check_valid=], [quality=], [firstface_index=], [flip_faces=])
 808// Description:
 809//   Takes a 2D input path, polygon or region and returns a path offset by the specified amount.  As with the built-in
 810//   offset() module, you can use `r` to specify rounded offset and `delta` to specify offset with
 811//   corners.  If you used `delta` you can set `chamfer` to true to get chamfers.
 812//   For paths and polygons positive offsets make the polygons larger.  For paths, 
 813//   positive offsets shift the path to the left, relative to the direction of the path.
 814//   .
 815//   If you use `delta` without chamfers, the path must not include any 180 degree turns, where the path
 816//   reverses direction.  Such reversals result in an offset with two parallel segments, so they cannot be
 817//   extended to an intersection point.  If you select chamfering the reversals are permitted and will result
 818//   in a single segment connecting the parallel segments.  With rounding, a semi-circle will connect the two offset segments.
 819//   Note also that repeated points are always illegal in the input; remove them first with {{deduplicate()}}.  
 820//   .
 821//   When offsets shrink the path, segments cross and become invalid.  By default `offset()` checks
 822//   for this situation.  To test validity the code checks that segments have distance larger than (r
 823//   or delta) from the input path.  This check takes O(N^2) time and may mistakenly eliminate
 824//   segments you wanted included in various situations, so you can disable it if you wish by setting
 825//   check_valid=false.  When segments are mistakenly removed, you may get the wrong offset output, or you may
 826//   get an error, depending on the effect of removing the segment.  
 827//   The erroneous removal of segments is more common when your input
 828//   contains very small segments and in this case can result in an invalid situation where the remaining
 829//   valid segments are parallel and cannot be connected to form an offset curve.  If this happens, you
 830//   will get an error message to this effect.  The only solutions are to either remove the small segments with {{deduplicate()}},
 831//   or if your path permits it, to set check_valid to false.  
 832//   .
 833//   Another situation that can arise with validity testing is that the test is not sufficiently thorough and some
 834//   segments persist that should be eliminated.  In this case, increase `quality` from its default of 1 to a value of 2 or 3.
 835//   This increases the number of samples on the segment that are checked, so it will increase run time.  In
 836//   some situations you may be able to decrease run time by setting quality to 0, which causes only
 837//   segment ends to be checked.  
 838//   .
 839//   When invalid segments are eliminated, the path length decreases, and multiple points on the input path map to the same point
 840//   on the offset path.  If you use chamfering or rounding, then
 841//   the chamfers and roundings can increase the length of the output path.  Hence points in the output may be 
 842//   difficult to associate with the input.  If you want to maintain alignment between the points you
 843//   can use the `same_length` option.  This option requires that you use `delta=` with `chamfer=false` to ensure
 844//   that no points are added.  with `same_length`, when points collapse to a single point in the offset, the output includes
 845//   that point repeated to preserve the correct length.  Generally repeated points will not appear in the offset output
 846//   unless you set `same_length` to true, but in some rare circumstances involving very short segments, it is possible for the
 847//   repeated points to occur in the output, even when `same_length=false`.  
 848//   .
 849//   Another way to obtain alignment information is to use the return_faces option, which can
 850//   provide alignment information for all offset parameters: it returns a face list which lists faces between
 851//   the original path and the offset path where the vertices are ordered with the original path
 852//   first, starting at `firstface_index` and the offset path vertices appearing afterwords.  The
 853//   direction of the faces can be flipped using `flip_faces`.  When you request faces the return
 854//   value is a list: [offset_path, face_list].
 855// Arguments:
 856//   path = the path to process.  A list of 2d points.
 857//   ---
 858//   r = offset radius.  Distance to offset.  Will round over corners.
 859//   delta = offset distance.  Distance to offset with pointed corners.
 860//   chamfer = chamfer corners when you specify `delta`.  Default: false
 861//   closed = if true path is treate as a polygon. Default: False.
 862//   check_valid = perform segment validity check.  Default: True.
 863//   quality = validity check quality parameter, a small integer.  Default: 1.
 864//   same_length = return a path with the same length as the input.  Only compatible with `delta=`.  Default: false
 865//   return_faces = return face list.  Default: False.
 866//   firstface_index = starting index for face list.  Default: 0.
 867//   flip_faces = flip face direction.  Default: false
 868// Example(2D,NoAxes): Offset the red star out by 10 units.
 869//   star = star(5, r=100, ir=30);
 870//   stroke(closed=true, star, width=3, color="red");
 871//   stroke(closed=true, width=3, offset(star, delta=10, closed=true));
 872// Example(2D,NoAxes):  Offset the star with chamfering
 873//   star = star(5, r=100, ir=30);
 874//   stroke(closed=true, star, width=3, color="red");
 875//   stroke(closed=true, width=3,
 876//          offset(star, delta=10, chamfer=true, closed=true));
 877// Example(2D,NoAxes):  Offset the star with rounding
 878//   star = star(5, r=100, ir=30);
 879//   stroke(closed=true, star, width=3, color="red");
 880//   stroke(closed=true, width=3,
 881//          offset(star, r=10, closed=true));
 882// Example(2D,NoAxes): Offset inward 
 883//   star = star(7, r=120, ir=50);
 884//   stroke(closed=true, width=3, star, color="red");
 885//   stroke(closed=true, width=3,
 886//          offset(star, delta=-15, closed=true));
 887// Example(2D,NoAxes): Inward offset with chamfers
 888//   star = star(7, r=120, ir=50);
 889//   stroke(closed=true, width=3, star, color="red");
 890//   stroke(closed=true, width=3,
 891//          offset(star, delta=-15, chamfer=true, closed=true));
 892// Example(2D,NoAxes):  Inward offset with rounding
 893//   star = star(7, r=120, ir=50);
 894//   stroke(closed=true, width=3, star, color="red");
 895//   stroke(closed=true, width=3,
 896//          offset(star, r=-15, closed=true, $fn=20));
 897// Example(2D): Open path.  The red path moves from left to right as shown by the arrow and the positive offset shifts to the left of the initial red path.
 898//   sinpath = 2*[for(theta=[-180:5:180]) [theta/4,45*sin(theta)]];
 899//   stroke(sinpath, width=2, color="red", endcap2="arrow2");
 900//   stroke(offset(sinpath, r=17.5),width=2);
 901// Example(2D,NoAxes): An open path in red with with its positive offset in yellow and its negative offset in blue. 
 902//   seg = [[0,0],[0,50]];
 903//   stroke(seg,color="red",endcap2="arrow2"); 
 904//   stroke(offset(seg,r=15,closed=false));
 905//   stroke(offset(seg,r=-15,closed=false),color="blue");
 906// Example(2D,NoAxes): Offsetting the same line segment closed=true.  On the left, we use delta with chamfer=false, in the middle, chamfer=true, and on the right, rounding with r=.  A "closed" path here means that the path backtracks over itself.  When this happens, a flat end occurs in the first case, an end with chamfered corners if chamfering is on, or a semicircular rounding in the rounded case.  
 907//   seg = [[0,0],[0,50]];
 908//   stroke(seg,color="red"); 
 909//   stroke([offset(seg,delta=15,closed=true)]);
 910//   right(45){
 911//     stroke(seg,color="red");
 912//     stroke([offset(seg,delta=15,chamfer=true,closed=true)]);
 913//   }
 914//   right(90){
 915//     stroke(seg,color="red");
 916//     stroke([offset(seg,r=15,closed=true)]);
 917//   }
 918// Example(2D,NoAxes): Cutting a notch out of a square with a path reversal
 919//   path = [[-10,-10],[-10,10],[0,10],[0,0],[0,10],[10,10],[10,-10]];
 920//   stroke([path],width=.25,color="red");
 921//   stroke([offset(path, r=-2,$fn=32,closed=true)],width=.25);
 922// Example(2D,NoAxes): A more complex example where the path turns back on itself several times.  
 923//   $fn=32;
 924//   path = [
 925//           [0,0], [5,5],
 926//           [10,0],[5,5],
 927//           [11,8],[5,5],
 928//           [5,10],[5,5],
 929//           [-1,4],[5,5]
 930//           ];
 931//   op=offset(path, r=1.5,closed=true);
 932//   stroke([path],width=.1,color="red");
 933//   stroke([op],width=.1);
 934// Example(2D,NoAxes):  With the default quality value, this case produces the wrong answer.  This happens because the offset edge corresponding to the long left edge (shown in green) is erroneously flagged as invalid.  If you use `r=` instead of `delta=` then this will fail with an error.  
 935//   test = [[0,0],[10,0],[10,7],[0,7], [-1,-3]];
 936//   polygon(offset(test,delta=-1.9, closed=true)); 
 937//   stroke([test],width=.1,color="red");
 938//   stroke(select(test,-2,-1), width=.1, color="green");
 939// Example(2D,NoAxes):  Using `quality=2` produces the correct result
 940//   test = [[0,0],[10,0],[10,7],[0,7], [-1,-3]];
 941//   polygon(offset(test,r=-1.9, closed=true, quality=2));   
 942//   stroke([test],width=.1,color="red");
 943// Example(2D,NoAxes): This case fails if `check_valid=true` when delta is large enough because segments are too close to the opposite side of the curve so they all get flagged as invalid and deleted from the output.  
 944//   star = star(5, r=22, ir=13);
 945//   stroke(star,width=.3,closed=true);                                                           
 946//   color("green")
 947//     stroke(offset(star, delta=-9, closed=true),width=.3,closed=true); // Works with check_valid=true (the default)
 948//   color("red")
 949//     stroke(offset(star, delta=-10, closed=true, check_valid=false),   // Fails if check_valid=true 
 950//            width=.3,closed=true); 
 951// Example(2D): But if you use rounding with offset then you need `check_valid=true` when `r` is big enough.  It works without the validity check as long as the offset shape retains a some of the straight edges at the star tip, but once the shape shrinks smaller than that, it fails.  There is no simple way to get a correct result for the case with `r=10`, because as in the previous example, it will fail if you turn on validity checks.  
 952//   star = star(5, r=22, ir=13);
 953//   color("green")
 954//     stroke(offset(star, r=-8, closed=true,check_valid=false), width=.1, closed=true);
 955//   color("red")
 956//     stroke(offset(star, r=-10, closed=true,check_valid=false), width=.1, closed=true);
 957// Example(2D,NoAxes): The extra triangles in this example show that the validity check cannot be skipped
 958//   ellipse = scale([20,4], p=circle(r=1,$fn=64));
 959//   stroke(ellipse, closed=true, width=0.3);
 960//   stroke(offset(ellipse, r=-3, check_valid=false, closed=true),
 961//          width=0.3, closed=true);
 962// Example(2D,NoAxes): The triangles are removed by the validity check
 963//   ellipse = scale([20,4], p=circle(r=1,$fn=64));
 964//   stroke(ellipse, closed=true, width=0.3);
 965//   stroke(offset(ellipse, r=-3, check_valid=true, closed=true),
 966//          width=0.3, closed=true);
 967// Example(2D,NoAxes): The region shown in red has the yellow offset region. 
 968//   rgn = difference(circle(d=100),
 969//                    union(square([20,40], center=true),
 970//                          square([40,20], center=true)));
 971//   stroke(rgn, width=1, color="red");
 972//   region(offset(rgn, r=-5));
 973// Example(2D,NoAxes): Using `same_length=true` to align the original curve to the offset.  Note that lots of points map to the corner at the top.
 974//   closed=false;
 975//   path = [for(angle=[0:5:180]) 10*[angle/100,2*sin(angle)]];
 976//   opath = offset(path, delta=-3,same_length=true,closed=closed);
 977//   stroke(path,closed=closed,width=.3);
 978//   stroke(opath,closed=closed,width=.3);
 979//   color("red") for(i=idx(path)) stroke([path[i],opath[i]],width=.3);
 980
 981function offset(
 982    path, r=undef, delta=undef, chamfer=false,
 983    closed=false, check_valid=true,
 984    quality=1, return_faces=false, firstface_index=0,
 985    flip_faces=false, same_length=false
 986) =
 987    assert(!(same_length && return_faces), "Cannot combine return_faces with same_length")
 988    is_region(path)?
 989        assert(!return_faces, "return_faces not supported for regions.")
 990        let(
 991            ofsregs = [for(R=region_parts(path))
 992                difference([for(i=idx(R)) offset(R[i], r=u_mul(i>0?-1:1,r), delta=u_mul(i>0?-1:1,delta),
 993                                      chamfer=chamfer, check_valid=check_valid, quality=quality,same_length=same_length,closed=true)])]
 994        )
 995        union(ofsregs)
 996    :
 997    let(rcount = num_defined([r,delta]))
 998    assert(rcount==1,"Must define exactly one of 'delta' and 'r'")
 999    assert(!same_length || (is_def(delta) && !chamfer), "Must specify delta, with chamfer=false, when same_length=true")
1000    assert(is_path(path), "Input must be a path or region")
1001    let(
1002        chamfer = is_def(r) ? false : chamfer,
1003        quality = max(0,round(quality)),
1004        flip_dir = closed && !is_polygon_clockwise(path)? -1 : 1,
1005        d = flip_dir * (is_def(r) ? r : delta)
1006    )
1007    d==0 && !return_faces ? path :
1008    let(
1009        shiftsegs = [for(i=[0:len(path)-2]) _shift_segment([path[i],path[i+1]], d),
1010                     if (closed) _shift_segment([last(path),path[0]],d)
1011                     else [path[0],path[1]]  // dummy segment, not used
1012                    ],
1013        // good segments are ones where no point on the segment is less than distance d from any point on the path
1014        good = check_valid ? _good_segments(path, abs(d), shiftsegs, closed, quality)
1015                           : repeat(true,len(shiftsegs)),
1016        goodsegs = bselect(shiftsegs, good),
1017        goodpath = bselect(path,good)
1018    )
1019    assert(len(goodsegs)-(!closed && select(good,-1)?1:0)>0,"Offset of path is degenerate")
1020    let(
1021        // Extend the shifted segments to their intersection points.  For open curves the endpoints
1022        // are simply the endpoints of the shifted segments.  If segments are parallel then the intersection
1023        // points will be undef
1024        sharpcorners = [for(i=[0:len(goodsegs)-1])
1025                             !closed && i==0 ? goodsegs[0][0]
1026                           : !closed && i==len(goodsegs)-1 ? goodsegs[len(goodsegs)-2][1]
1027                           : _segment_extension(select(goodsegs,i-1), select(goodsegs,i))],
1028
1029        // true if sharpcorner has two parallel segments that go in the same direction 
1030        cornercheck = [for(i=idx(goodsegs)) (!closed && (i==0 || i==len(goodsegs)-1))
1031                                          || is_def(sharpcorners[i])
1032                                          || approx(unit(deltas(select(goodsegs,i-1))[0]) * unit(deltas(goodsegs[i])[0]),-1)],
1033        dummyA = assert(len(sharpcorners)==2 || all(cornercheck),"Two consecutive valid offset segments are parallel but do not meet at their ends, maybe because path contains very short segments that were mistakenly flagged as invalid; unable to compute offset"),
1034        reversecheck = 
1035            !same_length 
1036              || !(is_def(delta) && !chamfer)            // Reversals only a problem in delta mode without chamfers
1037              || all_defined(sharpcorners),
1038        dummyB = assert(reversecheck, "Either validity check failed and removed a valid segment or the input 'path' contains a segment that reverses direction (180 deg turn).  Path reversals are not allowed when same_length is true because they increase path length."),
1039        // This is a Boolean array that indicates whether a corner is an outside or inside corner
1040        // For outside corners, the new corner is an extension (angle 0), for inside corners, it turns backward (angle 180)
1041        // If either side turns back it is an inside corner---must check both.
1042        // Outside corners can get rounded (if r is specified and there is space to round them)
1043        // We flag endpoints of open paths as inside corners so that we don't try to round
1044        outsidecorner =
1045            len(sharpcorners)==2 ? [closed,closed]
1046          : [for(i=idx(goodsegs))
1047                !closed && (i==0 || i==len(goodsegs)-1) ? false  // endpoints of open path never get rounded
1048              : is_undef(sharpcorners[i]) ? true
1049              : let(prevseg=select(goodsegs,i-1))
1050                (goodsegs[i][1]-goodsegs[i][0]) * (goodsegs[i][0]-sharpcorners[i]) > 0
1051                  && (prevseg[1]-prevseg[0]) * (sharpcorners[i]-prevseg[1]) > 0
1052            ],
1053        steps = is_def(delta) ? undef
1054              : [
1055                 for(i=[0:len(goodsegs)-1])  
1056                    r==0 ? 0
1057                  : !closed && (i==0 || i==len(goodsegs)-1) ? 0    // We don't round ends of open paths
1058                     // floor is important here to ensure we don't generate extra segments when nearly straight paths expand outward
1059                  : let(vang = vector_angle(select(goodsegs,i-1)[1]-goodpath[i],
1060                                            goodsegs[i][0]-goodpath[i]))
1061                    assert(!outsidecorner[i] || vang!=0,    // If outsidecorner[i] is true then vang>0 needed to give valid step count
1062                           "Offset computation failed, probably because validity check mistakenly removed a valid segment.  Increasing quality might fix this.")
1063                    1+floor(segs(r)*vang/360)
1064                ],
1065        // newcorners is a list where each entry is a list of the points that correspond to a single point in the sharpcorners 
1066        // list: newcorners[i] is the point list that replaces goodpath[i].  Without rounding or chamfering (or reversals),
1067        // this means each entry of newcorners is a singleton list.  But in the other cases, multiple points may appear at
1068        // a given position; newcorners later gets flattened to produce the final list, but the structure is needed to
1069        // establish point alignment for creating faces, or for duplicating points if same_length is true.  
1070        newcorners =
1071            [for(i=idx(goodsegs))
1072                 let(
1073                     basepts = [ select(goodsegs,i-1)[1], goodsegs[i][0] ]
1074                 )
1075                 is_def(sharpcorners[i]) &&
1076                   ((is_def(steps) && steps[i] <=1)  // Don't round if steps is smaller than 2
1077                     || !outsidecorner[i])           // Don't round inside corners
1078                ? [sharpcorners[i]]
1079                : chamfer ? _offset_chamfer(
1080                                  goodpath[i], [
1081                                      select(goodsegs,i-1)[1],
1082                                      sharpcorners[i],
1083                                      goodsegs[i][0]
1084                                  ], d
1085                              )
1086                : is_def(delta) ?
1087                      (
1088                         is_def(sharpcorners[i]) ? [sharpcorners[i]]
1089                       : let(normal = d*line_normal(basepts))
1090                         basepts + [normal,normal]
1091                      )
1092                : // rounded case
1093                  let(
1094                      class =_tri_class( [ each select(goodsegs,i-1), goodsegs[i][0]]),
1095                      cw = class==1,
1096                      ccw = class==-1
1097                  )
1098                  arc(cp=goodpath[i], cw=cw, ccw=ccw,
1099                      points=basepts,
1100                      n=steps[i]+3)
1101              ],
1102        pointcount = [for(entry=newcorners) len(entry)],
1103        edges = flatten(newcorners),
1104        faces = !return_faces? []
1105              : _makefaces(
1106                           flip_faces, firstface_index, good,
1107                           pointcount, closed
1108                          ),
1109        final_edges = same_length ? select(edges,
1110                                           [0,
1111                                            each list_head(cumsum([for(g=good) g?1:0]))
1112                                           ]
1113                                    )
1114                                  : edges
1115    ) return_faces? [edges,faces] : final_edges;
1116
1117
1118
1119/// Internal Function: _filter_region_parts()
1120///
1121/// splits region1 into subpaths where either it touches itself or crosses region2.  Classifies all of the
1122/// subpaths as described below and keeps the ones listed in keep1.  A similar process is performed for region2.
1123/// All of the kept subpaths are assembled into polygons and returned as a lst.
1124/// .
1125/// The four types of subpath from the region are defined relative to the second region:
1126///    "O" - the subpath is outside the second region
1127///    "I" - the subpath is in the second region's interior
1128///    "S" - the subpath is on the 2nd region's border and the two regions interiors are on the same side of the subpath
1129///    "U" - the subpath is on the 2nd region's border and the two regions meet at the subpath from opposite sides
1130/// You specify which type of subpaths to keep with a string of the desired types such as "OS".  
1131function _filter_region_parts(region1, region2, keep, eps=EPSILON) = 
1132    // We have to compute common vertices between paths in the region because
1133    // they can be places where the path must be cut, even though they aren't
1134    // found my the split_path function.  
1135    let(
1136        subpaths = split_region_at_region_crossings(region1,region2,eps=eps),
1137        regions=[force_region(region1),
1138                 force_region(region2)]
1139    )        
1140    _assemble_path_fragments(
1141        [for(i=[0:1])
1142           let(
1143               keepS = search("S",keep[i])!=[],
1144               keepU = search("U",keep[i])!=[],        
1145               keepoutside = search("O",keep[i]) !=[],
1146               keepinside = search("I",keep[i]) !=[],
1147               all_subpaths = flatten(subpaths[i])
1148           )
1149           for (subpath = all_subpaths)
1150               let(
1151                   midpt = mean([subpath[0], subpath[1]]),
1152                   rel = point_in_region(midpt,regions[1-i],eps=eps),
1153                   keepthis = rel<0 ? keepoutside
1154                            : rel>0 ? keepinside
1155                            : !(keepS || keepU) ? false
1156                            : let(
1157                                  sidept = midpt + 0.01*line_normal(subpath[0],subpath[1]),
1158                                  rel1 = point_in_region(sidept,regions[0],eps=eps)>0,
1159                                  rel2 = point_in_region(sidept,regions[1],eps=eps)>0
1160                              )
1161                              rel1==rel2 ? keepS : keepU
1162               )
1163               if (keepthis) subpath
1164        ]
1165    );
1166
1167
1168function _list_three(a,b,c) =
1169   is_undef(b) ? a : 
1170   [
1171     a,
1172     if (is_def(b)) b,
1173     if (is_def(c)) c
1174   ];
1175
1176
1177
1178// Function&Module: union()
1179// Synopsis: Performs a Boolean union operation.
1180// SynTags: Geom, Region
1181// Topics: Boolean Operations, Regions, Polygons, Shapes2D, Shapes3D
1182// See Also: difference(), intersection(), diff(), intersect(), exclusive_or()
1183// Usage:
1184//   union() CHILDREN;
1185//   region = union(regions);
1186//   region = union(REGION1,REGION2);
1187//   region = union(REGION1,REGION2,REGION3);
1188// Description:
1189//   When called as a function and given a list of regions or 2D polygons,
1190//   returns the union of all given regions and polygons.  Result is a single region.
1191//   When called as the built-in module, makes the union of the given children.
1192// Arguments:
1193//   regions = List of regions to union.
1194// Example(2D):
1195//   shape1 = move([-8,-8,0], p=circle(d=50));
1196//   shape2 = move([ 8, 8,0], p=circle(d=50));
1197//   color("green") region(union(shape1,shape2));
1198//   for (shape = [shape1,shape2]) color("red") stroke(shape, width=0.5, closed=true);
1199function union(regions=[],b=undef,c=undef,eps=EPSILON) =
1200    let(regions=_list_three(regions,b,c))
1201    len(regions)==0? [] :
1202    len(regions)==1? regions[0] :
1203    let(regions=[for (r=regions) is_path(r)? [r] : r])
1204    union([
1205           _filter_region_parts(regions[0],regions[1],["OS", "O"], eps=eps),           
1206           for (i=[2:1:len(regions)-1]) regions[i]
1207          ],
1208          eps=eps
1209    );
1210
1211
1212// Function&Module: difference()
1213// Synopsis: Performs a Boolean difference operation.
1214// SynTags: Geom, Region
1215// Topics: Boolean Operations, Regions, Polygons, Shapes2D, Shapes3D
1216// See Also: union(), intersection(), diff(), intersect(), exclusive_or()
1217// Usage:
1218//   difference() CHILDREN;
1219//   region = difference(regions);
1220//   region = difference(REGION1,REGION2);
1221//   region = difference(REGION1,REGION2,REGION3);
1222// Description:
1223//   When called as a function, and given a list of regions or 2D polygons, 
1224//   takes the first region or polygon and differences away all other regions/polygons from it.  The resulting
1225//   region is returned.
1226//   When called as the built-in module, makes the set difference of the given children.
1227// Arguments:
1228//   regions = List of regions or polygons to difference.
1229// Example(2D):
1230//   shape1 = move([-8,-8,0], p=circle(d=50));
1231//   shape2 = move([ 8, 8,0], p=circle(d=50));
1232//   for (shape = [shape1,shape2]) color("red") stroke(shape, width=0.5, closed=true);
1233//   color("green") region(difference(shape1,shape2));
1234function difference(regions=[],b=undef,c=undef,eps=EPSILON) =
1235     let(regions = _list_three(regions,b,c))
1236     len(regions)==0? []
1237   : len(regions)==1? regions[0]
1238   : regions[0]==[] ? []
1239   : let(regions=[for (r=regions) is_path(r)? [r] : r])
1240     difference([
1241                 _filter_region_parts(regions[0],regions[1],["OU", "I"], eps=eps),                
1242                 for (i=[2:1:len(regions)-1]) regions[i]
1243                ],
1244                eps=eps
1245     );
1246
1247
1248// Function&Module: intersection()
1249// Synopsis: Performs a Boolean intersection operation.
1250// SynTags: Geom, Region
1251// Topics: Boolean Operations, Regions, Polygons, Shapes2D, Shapes3D
1252// See Also: difference(), union(), diff(), intersect(), exclusive_or()
1253// Usage:
1254//   intersection() CHILDREN;
1255//   region = intersection(regions);
1256//   region = intersection(REGION1,REGION2);
1257//   region = intersection(REGION1,REGION2,REGION3);
1258// Description:
1259//   When called as a function, and given a list of regions or polygons returns the
1260//   intersection of all given regions.  Result is a single region.
1261//   When called as the built-in module, makes the intersection of all the given children.
1262// Arguments:
1263//   regions = List of regions to intersect.
1264// Example(2D):
1265//   shape1 = move([-8,-8,0], p=circle(d=50));
1266//   shape2 = move([ 8, 8,0], p=circle(d=50));
1267//   for (shape = [shape1,shape2]) color("red") stroke(shape, width=0.5, closed=true);
1268//   color("green") region(intersection(shape1,shape2));
1269function intersection(regions=[],b=undef,c=undef,eps=EPSILON) =
1270     let(regions = _list_three(regions,b,c))
1271     len(regions)==0 ? []
1272   : len(regions)==1? regions[0]
1273   : regions[0]==[] || regions[1]==[] ? []   
1274   : intersection([
1275                   _filter_region_parts(regions[0],regions[1],["IS","I"],eps=eps),                       
1276                   for (i=[2:1:len(regions)-1]) regions[i]
1277                  ],
1278                  eps=eps
1279     );
1280
1281
1282
1283// Function&Module: exclusive_or()
1284// Synopsis: Performs a Boolean exclusive-or operation.
1285// SynTags: Geom, Region
1286// Topics: Boolean Operations, Regions, Polygons, Shapes2D, Shapes3D
1287// See Also: union(), difference(), intersection(), diff(), intersect()
1288// Usage:
1289//   exclusive_or() CHILDREN;
1290//   region = exclusive_or(regions);
1291//   region = exclusive_or(REGION1,REGION2);
1292//   region = exclusive_or(REGION1,REGION2,REGION3);
1293// Description:
1294//   When called as a function and given a list of regions or 2D polygons, 
1295//   returns the exclusive_or of all given regions.  Result is a single region.
1296//   When called as a module, performs a Boolean exclusive-or of up to 10 children.  Note that when
1297//   the input regions cross each other the exclusive-or operator will produce shapes that
1298//   meet at corners (non-simple regions), which do not render in CGAL.  
1299// Arguments:
1300//   regions = List of regions or polygons to exclusive_or
1301// Example(2D): As Function.  A linear_sweep of this shape fails to render in CGAL.  
1302//   shape1 = move([-8,-8,0], p=circle(d=50));
1303//   shape2 = move([ 8, 8,0], p=circle(d=50));
1304//   for (shape = [shape1,shape2])
1305//       color("red") stroke(shape, width=0.5, closed=true);
1306//   color("green") region(exclusive_or(shape1,shape2));
1307// Example(2D): As Module.  A linear_extrude() of the resulting geometry fails to render in CGAL.  
1308//   exclusive_or() {
1309//       square(40,center=false);
1310//       circle(d=40);
1311//   }
1312function exclusive_or(regions=[],b=undef,c=undef,eps=EPSILON) =
1313     let(regions = _list_three(regions,b,c))
1314     len(regions)==0? []
1315   : len(regions)==1? force_region(regions[0])
1316   : regions[0]==[] ? exclusive_or(list_tail(regions))
1317   : regions[1]==[] ? exclusive_or(list_remove(regions,1))
1318   : exclusive_or([
1319                   _filter_region_parts(regions[0],regions[1],["IO","IO"],eps=eps),                  
1320                   for (i=[2:1:len(regions)-1]) regions[i]
1321                  ],
1322                  eps=eps
1323     );
1324
1325
1326module exclusive_or() {
1327    if ($children==1) {
1328        children();
1329    } else if ($children==2) {
1330        difference() {
1331            children(0);
1332            children(1);
1333        }
1334        difference() {
1335            children(1);
1336            children(0);
1337        }
1338    } else if ($children==3) {
1339        exclusive_or() {
1340            exclusive_or() {
1341                children(0);
1342                children(1);
1343            }
1344            children(2);
1345        }
1346    } else if ($children==4) {
1347        exclusive_or() {
1348            exclusive_or() {
1349                children(0);
1350                children(1);
1351            }
1352            exclusive_or() {
1353                children(2);
1354                children(3);
1355            }
1356        }
1357    } else if ($children==5) {
1358        exclusive_or() {
1359            exclusive_or() {
1360                children(0);
1361                children(1);
1362                children(2);
1363                children(3);
1364            }
1365            children(4);
1366        }
1367    } else if ($children==6) {
1368        exclusive_or() {
1369            exclusive_or() {
1370                children(0);
1371                children(1);
1372                children(2);
1373                children(3);
1374            }
1375            children(4);
1376            children(5);
1377        }
1378    } else if ($children==7) {
1379        exclusive_or() {
1380            exclusive_or() {
1381                children(0);
1382                children(1);
1383                children(2);
1384                children(3);
1385            }
1386            children(4);
1387            children(5);
1388            children(6);
1389        }
1390    } else if ($children==8) {
1391        exclusive_or() {
1392            exclusive_or() {
1393                children(0);
1394                children(1);
1395                children(2);
1396                children(3);
1397            }
1398            exclusive_or() {
1399                children(4);
1400                children(5);
1401                children(6);
1402                children(7);
1403            }
1404        }
1405    } else if ($children==9) {
1406        exclusive_or() {
1407            exclusive_or() {
1408                children(0);
1409                children(1);
1410                children(2);
1411                children(3);
1412            }
1413            exclusive_or() {
1414                children(4);
1415                children(5);
1416                children(6);
1417                children(7);
1418            }
1419            children(8);
1420        }
1421    } else if ($children==10) {
1422        exclusive_or() {
1423            exclusive_or() {
1424                children(0);
1425                children(1);
1426                children(2);
1427                children(3);
1428            }
1429            exclusive_or() {
1430                children(4);
1431                children(5);
1432                children(6);
1433                children(7);
1434            }
1435            children(8);
1436            children(9);
1437        }
1438    } else {
1439        assert($children<=10, "exclusive_or() can only handle up to 10 children.");
1440    }
1441}
1442
1443
1444
1445// Function&Module: hull_region()
1446// Synopsis: Compute convex hull of region or 2d path
1447// SynTags: Geom, Path
1448// Topics: Regions, Polygons, Shapes2D
1449// Usage:
1450//    path = hull_region(region);
1451//    hull_region(region);
1452// Description:
1453//   Given a path, or a region, compute the convex hull
1454//   and return it as a path.  This differs from {{hull()}} and {{hull2d_path()}} which
1455//   return an index list into the point list.  As a module invokes the native hull() on
1456//   the specified region.  
1457// Arguments:
1458//   region = region or path listing points to compute the hull from.  
1459// Example(2D, NoAxes):
1460//   data = [star(id=10,od=20,n=9),
1461//           right(30, star(id=12,od=25, n=7))];
1462//   stroke(data);
1463//   stroke([hull_region(data)],color="red");
1464function hull_region(region) =
1465  assert(is_path(region) || is_region(region))
1466  let(
1467      pts = is_region(region) ? flatten(region)
1468                              : region,
1469      order = hull2d_path(pts)
1470  )
1471  select(pts,order);
1472
1473module hull_region(region)
1474{
1475  hull()region(region);
1476}
1477
1478
1479// vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap